Annotation¶

Dendritic and macrophage spleenic subsets. Progenitor

Lewis et al., 2011

In [1]:
options(warn=-1)
In [2]:
library_load <- suppressMessages(
    
    list(
        
        # Seurat 
        library(Seurat), 
        
        # Data 
        library(tidyverse), 
        
        # Plot 
        library(ggplot2), 
        library(ComplexHeatmap), 
        library(ggplotify), 
        library(grid), 
        library(circlize), 
        
        # Python 
        library(reticulate)
        
    )
)
In [3]:
random_seed <- 42
set.seed(random_seed)
In [4]:
# Configure reticulate 
use_condaenv(condaenv='p.3.8.12-FD20200109SPLENO', conda="/nobackup/peer/fdeckert/miniconda3/bin/conda", required=NULL)
py_config()

python:         /nobackup/peer/fdeckert/miniconda3/envs/p.3.8.12-FD20200109SPLENO/bin/python
libpython:      /nobackup/peer/fdeckert/miniconda3/envs/p.3.8.12-FD20200109SPLENO/lib/libpython3.8.so
pythonhome:     /nobackup/peer/fdeckert/miniconda3/envs/p.3.8.12-FD20200109SPLENO:/nobackup/peer/fdeckert/miniconda3/envs/p.3.8.12-FD20200109SPLENO
version:        3.8.12 | packaged by conda-forge | (default, Jan 30 2022, 23:42:07)  [GCC 9.4.0]
numpy:          /nobackup/peer/fdeckert/miniconda3/envs/p.3.8.12-FD20200109SPLENO/lib/python3.8/site-packages/numpy
numpy_version:  1.21.6

NOTE: Python version was forced by use_python function
In [5]:
ht_opt$message=FALSE # ComplexHeatmap
In [6]:
# Set working directory to project root
setwd("/research/peer/fdeckert/FD20200109SPLENO")
In [7]:
# Source files
source("plotting_global.R")
source("bin/cell_type.R")
source("bin/seurat_qc.R")
source("bin/pbDEA.R")
source("bin/raceid.R")
source("bin/gene_modules.R")
In [8]:
# Plotting Theme
ggplot2::theme_set(theme_global_set()) # From project global source()

Parameter settings¶

In [9]:
so_file <- "data/object/sct/so_sct_int_hvg8000.rds"
so_pp_file <- "data/object/pp.rds"
h5ad_pp_file <- "data/object/pp.h5ad"

Import Seurat object¶

In [10]:
so <- readRDS(so_file)

Seurat dimensional reduction and clustering¶

In [11]:
DefaultAssay(so) <- "integrated"

# Cluster all cells 
so <- FindNeighbors(so, dims=1:10, k.param=20, verbose=FALSE)
so <- FindClusters(so, verbose=FALSE, resolution=1, algorithm=1, group.singletons=TRUE)
In [12]:
# Umap dimensional reduction
so <- RunUMAP(so, dims=1:15, n.neighbors=100, min.dist=1, spread=1, verbose=FALSE, umap.method="umap-learn")
In [13]:
options(repr.plot.width=30, repr.plot.height=10)

dplot_1 <- dplot(so, reduction="umap", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so, reduction="umap", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_3 <- dplot(so, reduction="umap", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot_1 <- fplot(so, reduction="umap", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot_2 <- fplot(so, reduction="umap", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")
fplot_3 <- fplot(so, reduction="umap", features="Ccr2", slot="data") + ggtitle("Ccr2") + scale_color_viridis(option="G")
fplot_4 <- fplot(so, reduction="umap", features="Ccr7", slot="data") + ggtitle("Ccr7") + scale_color_viridis(option="G")
fplot_5 <- fplot(so, reduction="umap", features="Ly6c2", slot="data") + ggtitle("Ly6c2") + scale_color_viridis(option="G")

dplot_1 + dplot_2 + dplot_3 + fplot_1 + fplot_2 + fplot_3 + fplot_4 + fplot_5 + plot_layout(nrow=2, ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
In [14]:
so@meta.data %>% dplyr::group_by(seurat_clusters, treatment, sample_rep) %>% 
    dplyr::summarise(n=n()) %>% data.frame() %>% 
    tidyr::spread(seurat_clusters, n) %>% 
    kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()

`summarise()` has grouped output by 'seurat_clusters', 'treatment'. You can
override using the `.groups` argument.
treatment sample_rep 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
CpG Rep1 705 590 154 476 412 167 383 337 258 254 157 154 102 56 133 117 100 70 101 74 13 NA
CpG Rep2 930 654 160 648 108 348 106 574 226 450 343 314 296 150 293 98 66 161 38 20 43 36
NaCl Rep1 170 161 327 72 372 106 437 23 321 49 121 77 59 133 70 65 121 42 106 109 15 NA
NaCl Rep2 257 251 685 90 224 477 91 52 138 167 252 303 314 426 117 73 60 61 31 33 25 4
In [15]:
# Rmove cluster wich have no evidence in replicates 
so <- subset(so, subset=seurat_clusters!="21")
In [16]:
DefaultAssay(so) <- "RNA"

RaceID¶

In [17]:
raceid <- raceid_pp(so=so, suffix="", compute=FALSE)
so <- raceid_to_seurat(so, raceid)
In [18]:
options(repr.plot.width=6*6, repr.plot.height=6)

dplot_1 <- dplot(so, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
dplot_5 <- dplot(so, reduction="umap_varid", group_by="label_fine_haemosphere", alpha=0.5) + scale_color_manual(values=color$label_fine_haemosphere)
fplot_1 <- fplot(so, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")

dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(nrow=1, ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
In [19]:
so@meta.data %>% dplyr::group_by(varid_clusters, treatment, sample_rep) %>% 
    dplyr::summarise(n=n()) %>% data.frame() %>% 
    tidyr::spread(varid_clusters, n) %>% 
    kableExtra::kable("html") %>% as.character() %>% IRdisplay::display_html()

`summarise()` has grouped output by 'varid_clusters', 'treatment'. You can
override using the `.groups` argument.
treatment sample_rep 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
CpG Rep1 904 270 682 155 565 519 464 326 157 195 186 89 67 45 59 41 31 42 12 3 1
CpG Rep2 1294 653 984 207 266 262 664 447 340 183 277 64 83 94 29 88 62 10 13 5 1
NaCl Rep1 153 204 116 404 551 458 155 49 118 305 73 125 37 41 65 10 26 27 11 19 9
NaCl Rep2 445 935 235 893 171 306 209 107 246 101 120 61 47 49 28 41 37 41 26 16 13

Progenitor cell annotation¶

Progenitor gene modules¶

In [20]:
seurat_clusters_prog <- c("17", "15", "2", "13", "5", "12", "11", "9", "7", "3", "0", "1")
so_prog <- subset(so, subset=seurat_clusters %in% seurat_clusters_prog)
In [21]:
so_prog <- AddModuleScore(so_prog, list(genes_hsc), assay="RNA", name="msHSC")
so_prog <- AddModuleScore(so_prog, list(genes_ly), assay="RNA", name="msLy")
so_prog <- AddModuleScore(so_prog, list(genes_meg), assay="RNA", name="msMeg")
so_prog <- AddModuleScore(so_prog, list(genes_ery), assay="RNA", name="msEry")
so_prog <- AddModuleScore(so_prog, list(genes_mo), assay="RNA", name="msMo")
so_prog <- AddModuleScore(so_prog, list(genes_neu), assay="RNA", name="msNeu")
so_prog <- AddModuleScore(so_prog, list(genes_eo), assay="RNA", name="msEo")
so_prog <- AddModuleScore(so_prog, list(genes_baso), assay="RNA", name="msBaso")
so_prog <- AddModuleScore(so_prog, list(genes_mast), assay="RNA", name="msMast")

RaceID (Progenitors)¶

In [22]:
raceid_prog <- raceid_pp(so=so_prog, suffix="_prog", compute=FALSE)
so_prog <- raceid_to_seurat(so_prog, raceid_prog)
In [23]:
options(repr.plot.width=6*6, repr.plot.height=6)

dplot_1 <- dplot(so_prog, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_prog, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_prog, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_prog, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
dplot_5 <- dplot(so_prog, reduction="umap_varid", group_by="label_fine_haemosphere", alpha=0.5) + scale_color_manual(values=color$label_fine_haemosphere)
fplot_1 <- fplot(so_prog, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")

dplot_1 + dplot_2 + dplot_3 + dplot_4 + dplot_5 + fplot_1 + plot_layout(nrow=1, ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Marker gene heat map (Progenitors)¶

In [24]:
options(repr.plot.width=7.5, repr.plot.height=15)

# Get expression matrix 
mat <- AverageExpression(so_prog, assay="RNA", slot="data", features=genes_marker$genes, group.by="seurat_clusters")[["RNA"]]

mat <- mat[, seurat_clusters_prog]
mat <- t(scale(t(mat)))
breaks <- seq(-max(abs(mat)), +max(abs(mat)), by=0.1)

hm_1 <- ComplexHeatmap::Heatmap(
    
    matrix=mat,    
    name="z-score", 
    column_title="Progenitor marker expression",
    col=colorRamp2(breaks, mako(length(breaks))), 
    
    column_title_gp=gpar(fontsize=18, fontface="bold"), 
    column_names_gp =grid::gpar(fontsize=16), 
    row_title_gp=gpar(fontsize=16, fontface="bold"),
    row_names_gp=grid::gpar(fontsize=16), 
    
    row_split=genes_marker$cell_type,
    
    cluster_rows=TRUE, 
    cluster_columns=FALSE,
    show_row_names=TRUE,
    show_column_names=TRUE, 
    
    row_dend_width=unit(1.5, "cm"), 
    width=ncol(mat)*unit(5, "mm"), 
    height=nrow(mat)*unit(8, "mm"), 
    rect_gp=gpar(col="white", lwd=2), 
    
    heatmap_legend_param=list(title_gp=gpar(fontsize=16, fontface="bold"), labels_gp=gpar(fontsize=16))

) %>% as.ggplot()

hm_1

Annotate erythroblasts¶

In [25]:
seurat_clusters_eb <- c("13", "5", "12", "11", "9", "7", "3", "0", "1")
so_eb <- subset(so, subset=seurat_clusters %in% seurat_clusters_eb)
In [26]:
cell_type_eb <- data.frame(
    
    ident=c(
        
        "13", 
        "5", 
        "12", 
        "11", 
        "9", 
        "7", 
        "3", 
        "0", 
        "1"
        
    ), 
    
    cell_type_fine_detail=c(
    
        "ProEB (1)", 
        "ProEB (2)", 
        "ProEB (3)", 
        "ProEB (4)", 
        "EB (1)", 
        "EB (2)", 
        "EB (3)", 
        "EB (4)", 
        "EB (5)"
    
    ), 
    
    cell_type_fine=c(
        
        "ProEB (1)", 
        "ProEB (2)", 
        "ProEB (3)", 
        "ProEB (4)", 
        "EB (1)", 
        "EB (2)", 
        "EB (3)", 
        "EB (4)", 
        "EB (5)"
        
    ), 
    
    cell_type_main=c(
    
        "ProEB", 
        "ProEB",
        "ProEB",
        "ProEB",
        "EB", 
        "EB",
        "EB",
        "EB", 
        "EB"
        
    )
    
)
In [27]:
cell_type_eb <- dplyr::left_join(dplyr::select(so_eb@meta.data, seurat_clusters, cell_id), cell_type_eb, by=c("seurat_clusters"="ident")) %>% 
    column_to_rownames("cell_id")  %>% 
    dplyr::mutate(varid_clusters=paste0("varid_clusters_eb_NA")) %>% 
    dplyr::mutate(ident=paste0("seurat_clusters_eb_", seurat_clusters))

RaceID (MPP)¶

In [28]:
seurat_clusters_mpp <- c("17", "15", "2")
so_mpp <- subset(so_prog, subset=seurat_clusters %in% seurat_clusters_mpp)
In [29]:
raceid_mpp <- raceid_pp(so=so_mpp, suffix="_mpp", compute=FALSE)
so_mpp <- raceid_to_seurat(so_mpp, raceid_mpp)
In [30]:
options(repr.plot.width=6*6, repr.plot.height=3*6)

dplot_1 <- dplot(so_mpp, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_mpp, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_mpp, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_mpp, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot1 <- fplot(so_mpp, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot2 <- fplot(so_mpp, reduction="umap_varid", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")
fplot3 <- fplot(so_mpp, reduction="umap_varid", features="msHSC1") + ggtitle("msHSC1") + scale_color_viridis(option="G")
fplot4 <- fplot(so_mpp, reduction="umap_varid", features="msLy1") + ggtitle("msLy1") + scale_color_viridis(option="G")
fplot5 <- fplot(so_mpp, reduction="umap_varid", features="msMeg1") + ggtitle("msMeg1") + scale_color_viridis(option="G")
fplot6 <- fplot(so_mpp, reduction="umap_varid", features="msEry1") + ggtitle("msEry1") + scale_color_viridis(option="G")
fplot7 <- fplot(so_mpp, reduction="umap_varid", features="msMo1") + ggtitle("msMo1") + scale_color_viridis(option="G")
fplot8 <- fplot(so_mpp, reduction="umap_varid", features="msNeu1") + ggtitle("msNeu1") + scale_color_viridis(option="G")
fplot9 <- fplot(so_mpp, reduction="umap_varid", features="msEo1") + ggtitle("msEo1") + scale_color_viridis(option="G")
fplot10 <- fplot(so_mpp, reduction="umap_varid", features="msBaso1") + ggtitle("msBaso1") + scale_color_viridis(option="G")
fplot11 <- fplot(so_mpp, reduction="umap_varid", features="msMast1") + ggtitle("msMast1") + scale_color_viridis(option="G")

dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot1 + fplot2 + fplot3 + fplot4 + fplot5 + fplot6 + fplot7 + fplot8 + fplot9 + fplot10 + fplot11 + plot_layout(ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
In [31]:
options(repr.plot.width=7.5, repr.plot.height=15)

source("bin/gene_modules.R")

mat <- AverageExpression(so_mpp, assay="RNA", slot="data", features=genes_marker$genes, group.by="varid_clusters")[["RNA"]]

mat <- t(scale(t(mat)))
mat <- na.omit(mat) 

cell_type_mpp <- genes_marker$cell_type[genes_marker$genes %in% rownames(mat)]

breaks <- seq(-max(abs(mat)), +max(abs(mat)), by=0.1)

hm_1 <- ComplexHeatmap::Heatmap(
    
    matrix=mat,    
    name="z-score", 
    column_title="Progenitor marker expression",
    col=colorRamp2(breaks, mako(length(breaks))), 
    
    column_title_gp=gpar(fontsize=18, fontface="bold"), 
    column_names_gp =grid::gpar(fontsize=16), 
    row_title_gp=gpar(fontsize=16, fontface="bold"),
    row_names_gp=grid::gpar(fontsize=16), 
    
    row_split=cell_type_mpp,
    cluster_rows=TRUE, 
    cluster_columns=FALSE,
    show_row_names=TRUE,
    show_column_names=TRUE, 
    
    row_dend_width=unit(1.5, "cm"), 
    width=ncol(mat)*unit(5, "mm"), 
    height=nrow(mat)*unit(8, "mm"), 
    
    rect_gp=gpar(col="white", lwd=2), 
    heatmap_legend_param=list(title_gp=gpar(fontsize=16, fontface="bold"), labels_gp=gpar(fontsize=16))

) %>% as.ggplot()

hm_1

Annotate MPP¶

In [32]:
cell_type_mpp <- data.frame(
    
    ident=c(
        
        1, 
        2, 
        3, 
        4, 
        5, 
        6, 
        7, 
        8, 
        9, 
        10
        
    ), 
    
    cell_type_fine_detail=c(
    
        "MEP (2)", 
        "MEP (4)", 
        "MEP (3)", 
        "GMP", 
        "MegP", 
        "BasoP", 
        "MLP", 
        "MDP", 
        "MEP (1)", 
        "MastP"
    
    ), 
    
    cell_type_fine=c(
        
        "MEP (2)", 
        "MEP (4)", 
        "MEP (3)", 
        "GMP", 
        "MegP", 
        "BasoP", 
        "MLP", 
        "MDP", 
        "MEP (1)", 
        "MastP"

    ), 
    
    cell_type_main=c(
        
        "MEP", 
        "MEP", 
        "MEP", 
        "GMP", 
        "MegP", 
        "GMP", 
        "MLP", 
        "MDP", 
        "MEP", 
        "GMP"
    
    )

)
In [33]:
cell_type_mpp <- dplyr::left_join(dplyr::select(so_mpp@meta.data, seurat_clusters, varid_clusters, cell_id), cell_type_mpp, by=c("varid_clusters"="ident")) %>% 
    column_to_rownames("cell_id") %>% 
    dplyr::mutate(varid_clusters=paste0("varid_clusters_mpp_", varid_clusters)) %>% 
    dplyr::mutate(ident=varid_clusters)

Myeloid cell annotation¶

In [34]:
seurat_clusters_m <- c("10", "6", "18", "20", "19", "4", "8", "14", "16")
so_m <- subset(so, subset=seurat_clusters %in% seurat_clusters_m)

RaceID (Myeloid)¶

In [35]:
raceid_m <- raceid_pp(so=so_m, suffix="_m", compute=FALSE)
so_m <- raceid_to_seurat(so_m, raceid_m)
In [36]:
options(repr.plot.width=6*6, repr.plot.height=6)

dplot_1 <- dplot(so_m, reduction="umap_varid", group_by="seurat_clusters", alpha=0.5)
dplot_2 <- dplot(so_m, reduction="umap_varid", group_by="varid_clusters", alpha=0.5)
dplot_3 <- dplot(so_m, reduction="umap_varid", group_by="treatment", alpha=0.5) + scale_color_manual(values=color$treatment)
dplot_4 <- dplot(so_m, reduction="umap_varid", group_by="cc_phase_class", alpha=0.5) + scale_color_manual(values=color$cc_phase_class)
fplot_1 <- fplot(so_m, reduction="umap_varid", features="pRb_RNA") + ggtitle("Percentage Rb") + scale_color_viridis(option="G")
fplot_2 <- fplot(so_m, reduction="umap_varid", features="pMt_RNA") + ggtitle("Percentage Mt") + scale_color_viridis(option="G")

dplot_1 + dplot_2 + dplot_3 + dplot_4 + fplot_1 + fplot_2 + plot_layout(nrow=1, ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
In [37]:
options(repr.plot.width=6*6, repr.plot.height=6)

fplot_1 <- fplot(so_m, reduction="umap_varid", features="Ccr7", slot="data") + ggtitle("Ccr7 (TCZ)") + scale_color_viridis(option="G")
fplot_2 <- fplot(so_m, reduction="umap_varid", features="Cxcr5", slot="data") + ggtitle("Cxcr5 (BCZ)") + scale_color_viridis(option="G")
fplot_3 <- fplot(so_m, reduction="umap_varid", features="Cxcr4", slot="data") + ggtitle("Cxcr4 (RP)") + scale_color_viridis(option="G")
fplot_4 <- fplot(so_m, reduction="umap_varid", features="Ccr2", slot="data") + ggtitle("Ccr2 (Blood)") + scale_color_viridis(option="G")

fplot_1 + fplot_2 + fplot_3 + fplot_4 + plot_spacer() + plot_spacer() + plot_layout(nrow=1, ncol=6)

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.

Set varid_clusters to re-use the DEA results¶

In [38]:
so_m$varid_clusters <- paste0("varid_clusters_m_", so_m$varid_clusters)

Canonical marker genes (Myeloid)¶

In [39]:
features_m <- c(
    
    "Adgre1", # F4/20
    "Itgax", # Cd11c
    "Itgam", # Cd11b
    "Ly6c2", # Ly6c
    "Csf1r", # Cd115
    "Cd8a", # Cd8 
    "Cd4", # Cd4
    "Napsa", # Napsa
    "Lsp1" # Lsp1
    
)
In [40]:
options(repr.plot.width=3*(2.5+length(features_m)/2.5), repr.plot.height=5)

dp_1 <- dp_feature(so_m, features_m, group_by="varid_clusters", scale=FALSE, title="Myeloid marker", range_max=8) + theme_global_set(1)
dp_2 <- dp_feature(so_m, features_m, group_by="varid_clusters", scale=TRUE, title="Myeloid marker", range_max=8) + theme_global_set(1)

dp_1 + dp_2

`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
In [41]:
options(repr.plot.width=20, repr.plot.height=20)

results <- voomlmfit_marker(so_m, grouping_var="varid_clusters", pseudobatch_var="sample_group") 
hm_marker(results, so_m, grouping_var="varid_clusters", width=0.05, height=1)

 [1] "varid_clusters_m_1"  "varid_clusters_m_7"  "varid_clusters_m_6" 
 [4] "varid_clusters_m_3"  "varid_clusters_m_8"  "varid_clusters_m_5" 
 [7] "varid_clusters_m_4"  "varid_clusters_m_2"  "varid_clusters_m_10"
[10] "varid_clusters_m_11" "varid_clusters_m_9" 
 [1] "varid_clusters_m_1"  "varid_clusters_m_7"  "varid_clusters_m_6" 
 [4] "varid_clusters_m_3"  "varid_clusters_m_8"  "varid_clusters_m_5" 
 [7] "varid_clusters_m_4"  "varid_clusters_m_2"  "varid_clusters_m_10"
[10] "varid_clusters_m_11" "varid_clusters_m_9"

Macrophage marker genes¶

In [42]:
features_mac <- c(
    
    "Spic", # Spic | RPM marker
    "Vcam1", # Vcam1 | RPM marker
    "Slc40a1", # Slc40a1 | RPM marker | Iron transport
    "Siglec1", # Cd169 | MMM marker. Binds vaious cell type. 
    "Nr1h3", # Lxr alpha | MMMM marker. A transcription factor activated by oxysterols
    "Marco", # Marco | MZM marker
    "Cd209b", # SIGNR1 | MZM marker 
    "Gas6", # Gas6 | TBM
    "Mfge8"

)

split_mac <- c(

    "RPM", 
    "RPM", 
    "RPM",
    "MMM",
    "MMM", 
    "MZM", 
    "MZM", 
    "TBM", 
    "TBM"
    
)
In [43]:
options(repr.plot.width=3*(2.5+length(features_mac)/2.5), repr.plot.height=5)

dp_1 <- dp_feature(so_m, features_mac, split_mac, group_by="varid_clusters", scale=FALSE, title="Macrophage marker", range_max=8) + theme_global_set(1)
dp_2 <- dp_feature(so_m, features_mac, split_mac, group_by="varid_clusters", scale=TRUE, title="Macrophage marker", range_max=8) + theme_global_set(1)

dp_1 + dp_2

`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
In [44]:
options(repr.plot.width=10, repr.plot.height=12)

results <- voomlmfit_marker(subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_10", "varid_clusters_m_6")), grouping_var="varid_clusters", pseudobatch_var="sample_group") 
hm_marker(results, subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_10", "varid_clusters_m_6")), grouping_var="varid_clusters")

[1] "varid_clusters_m_6"  "varid_clusters_m_10"
[1] "varid_clusters_m_6"  "varid_clusters_m_10"

Dendritic cells marker genes¶

In [45]:
features_dc <- c(
    
    "Zbtb46", # Zbtb46 | Common DC progenitor (CDP) | TF
    "Xcr1", # Xcr1 | cDC1 of which most express CD8aa
    "Ly75", # Dec205 (Lectin receptor) | cDC1 subset in the WP | After immunization cDC1 colocalize with and preferentially activate Cd8+ T cells in the WP
    "Itgae", "Cd207", # Cd103 Langerin | cDC1 in the MZ and RP at stready-state
    "Sirpa", "Gpr183", "Mbtps1", # Sirpa Ebi2 (recognizes oxysterol ligands) S1P | cDC2 in the bridging channel (BC) at steady-state condition and also express Cd11b. Ebi2 and S1p maintain BC localization. 
    "Esam", "Notch2", "Rbpj", "Clec4a4", "Irf4", "Klf4", # Esam (endothelial cell adhesion moleculoe) Notch2 Rbpj Dcir2 | cDC2 Esam hi subset | Notch2 and Rbpj signaling is required for Esam expression. The subset also express Cd11b, Cd4 and Dcir2
    "Cx3cr1", "Il12a", "Il12b", "Tnf", # Cx3cr1 Il12 Il12 Infalpha | cDC2 Esam low subset | Independent of Notch2 and Irf4 signaling and producing innflammatory cytokines such as TNFa and Il12. Also positve for Cd11b, Dcir2 but less cd4 or double negative. 
    "Cd4", 
    "Cd8a", 
    "Ccr7", 
    "Cxcr5", 
    "Cxcr4", 
    "Ccr2", 
    "Pcna", 
    "Uhrf1"
    
)

split_dc <- c(

    "CDP", 
    rep("cDC1", 4), 
    rep("cDC2", 13), 
    rep("CD", 2), 
    rep("Chemokines", 4),
    rep("CC", 2)
    
)
In [46]:
options(repr.plot.width=2*(2.5+length(features_dc)/2.5), repr.plot.height=5)

dp_1 <- dp_feature(so_m, features_dc, split_dc, group_by="varid_clusters", scale=FALSE, range_max=8) + theme_global_set(1)
dp_2 <- dp_feature(so_m, features_dc, split_dc, group_by="varid_clusters", scale=TRUE, range_max=8) + theme_global_set(1)

dp_1 + dp_2

`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
In [47]:
options(repr.plot.width=10, repr.plot.height=12)

results <- voomlmfit_marker(subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_4", "varid_clusters_m_5", "varid_clusters_m_2", "varid_clusters_m_8", "varid_clusters_m_9")), grouping_var="varid_clusters", pseudobatch_var="sample_group") 
hm_marker(results, subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_4", "varid_clusters_m_5", "varid_clusters_m_2", "varid_clusters_m_8", "varid_clusters_m_9")), grouping_var="varid_clusters")

[1] "varid_clusters_m_8" "varid_clusters_m_5" "varid_clusters_m_4"
[4] "varid_clusters_m_2" "varid_clusters_m_9"
[1] "varid_clusters_m_8" "varid_clusters_m_5" "varid_clusters_m_4"
[4] "varid_clusters_m_2" "varid_clusters_m_9"

Monocyte marker genes¶

In [48]:
features_mo <- c(
    
    "Ly6c2", "Csf1r", "Spi1", "Spn", # Ly6c Ly6c Csf1r Pu1 Cd43 | Canonical marker genes | Two major monocyte populations in the blood which are Ly6c lo and Ly6c high. Cd43 is low in classical and high in int and non-classical. 
    "Jun", "Fos", "Irf8", "Klf4", "Vcan", "Cd163", "Cd63", "S100a8", # Ccr2 | Classical monocyte | Recruitment to blood by Ccr2
    "Treml4", "Cx3cr1", "Nr4a1", "Klf2", "Fcgr3", "Ifitm1", "Ifitm2", "Ifitm3", "Cdkn1c", "Mtss1", # Treml4 Cx3cr1 | Non-classical monocyte | Migrate in response to Cx3cr1 
    "H2-Ab1", "Cd74",  "Ccr5", "Cxcl9", "Cxcl10", 
    "Cd4", 
    "Cd8a", 
    "Ccr7", 
    "Cxcr5", 
    "Cxcr4", 
    "Ccr2", 
    "Pcna", 
    "Uhrf1"
)

split_mo <- c(

    rep("Canonical", 4), 
    rep("cMo", 8), 
    rep("ncMo", 10),
    rep("moDC", 5), 
    rep("CD", 2), 
    rep("Chemokines", 4),
    rep("CC", 2)
)
In [49]:
options(repr.plot.width=2*(2.5+length(features_mo)/2.5), repr.plot.height=5)

dp_1 <- dp_feature(so_m, features_mo, split_mo, group_by="varid_clusters", scale=FALSE, title="Monocyte marker", range_max=8) + theme_global_set(1)
dp_2 <- dp_feature(so_m, features_mo, split_mo, group_by="varid_clusters", scale=TRUE, title="Monocyte marker", range_max=8) + theme_global_set(1)

dp_1 + dp_2

`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
`summarise()` has grouped output by 'varid_clusters'. You can override using
the `.groups` argument.
Scale for fill is already present.
Adding another scale for fill, which will replace the existing scale.
In [50]:
options(repr.plot.width=10, repr.plot.height=12)

results <- voomlmfit_marker(subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_3", "varid_clusters_m_1", "varid_clusters_m_11", "varid_clusters_m_7")), grouping_var="varid_clusters", pseudobatch_var="sample_group") 
hm_marker(results, subset(so_m, subset=varid_clusters %in% c("varid_clusters_m_3", "varid_clusters_m_1", "varid_clusters_m_11", "varid_clusters_m_7")), grouping_var="varid_clusters")

[1] "varid_clusters_m_1"  "varid_clusters_m_7"  "varid_clusters_m_3" 
[4] "varid_clusters_m_11"
[1] "varid_clusters_m_1"  "varid_clusters_m_7"  "varid_clusters_m_3" 
[4] "varid_clusters_m_11"
In [51]:
cell_type_m <- data.frame(
    
    ident=c(
        
        1, 
        2, 
        3, 
        4, 
        5, 
        6, 
        7, 
        8, 
        9, 
        10, 
        11
        
    ), 
    
    cell_type_fine_detail=c(
        
        "ncMo Cd4- (1)", 
        "cDC2 (1)", 
        "cMo Ly6c lo (2)",
        "cDC1 Cd8+ prolif. (1)", 
        "cDC1 Cd8+ (2)", 
        "RPM", 
        "ncMo Cd4+ (2)", 
        "cDC2 Ccr7+ (3)",   
        "cDC2 prolif. (2)",  
        "PreRPM", 
        "cMo Ly6c hi (1)" 

    ),
    
    cell_type_fine=c(
        
        "ncMo (1)", 
        "cDC2 (1)", 
        "cMo (2)", 
        "cDC1 (1)",  
        "cDC1 (2)", 
        "RPM", 
        "ncMo (2)", 
        "cDC2 (3)",  
        "cDC2 (2)", 
        "PreRPM", 
        "cMo (1)" 

    ), 
    
    cell_type_main=c(
        
        "ncMo", 
        "cDC2", 
        "cMo", 
        "cDC1", 
        "cDC1", 
        "RPM", 
        "ncMo", 
        "cDC2", 
        "cDC2", 
        "PreRPM", 
        "cMo"
    
    )

)
In [52]:
cell_type_m <- dplyr::mutate(cell_type_m, ident=paste0("varid_clusters_m_", ident)) %>% 
    dplyr::left_join(dplyr::select(so_m@meta.data, seurat_clusters, varid_clusters, cell_id), ., by=c("varid_clusters"="ident")) %>% 
    column_to_rownames("cell_id") %>% 
    dplyr::mutate(ident=varid_clusters)

Combine annotation results¶

In [53]:
cell_type <- rbind(cell_type_eb, cell_type_mpp, cell_type_m) %>% dplyr::select(cell_type_main, cell_type_fine, cell_type_fine_detail, varid_clusters, ident)
so <- AddMetaData(so, cell_type)
In [54]:
unique(dplyr::select(cell_type, cell_type_main, cell_type_fine, cell_type_fine_detail, ident)) %>% write.csv("test.csv")
In [55]:
options(repr.plot.width=2*10, repr.plot.height=6)

# Source files
dplot_1 <- dplot(so, reduction="umap", group_by="cell_type_fine_detail", alpha=1, pt_size=1) + 
    scale_color_manual(values=color[["cell_type_fine_detail"]]) + 
    theme(
        legend.title=element_blank(), 
        plot.title=element_blank()
    )

dplot_2 <- dplot(so, reduction="umap", group_by="cell_type_main", alpha=1, pt_size=1) + 
    scale_color_manual(values=color[["cell_type_main"]]) + 
    theme(
        legend.title=element_blank(), 
        plot.title=element_blank()
    )

dplot_1 + dplot_2

Save results¶

In [56]:
# saveRDS(so, so_pp_file)
In [57]:
# write.csv(so@meta.data, "data/object/components/meta.csv")
# write.csv(rownames(so), "data/object/genes.csv", row.names=FALSE)
# write.csv(colnames(so), "data/object/cells.csv", row.names=FALSE)
# write.csv(so@reductions$umap@cell.embeddings, "data/object/components/umap.csv")
In [58]:
# # Store data as h5ad 
# adata <- import("anndata", as="adata", convert=FALSE)
# pd <- import("pandas", as="pd", convert=FALSE)
# np <- import("numpy", as="np", convert=FALSE)
    
# # Transform dgCMatrix to 
# X <- GetAssayData(so, assay="RNA", slot="counts") %>% as.matrix() %>% t()
# X <- np$array(X, dtype=np$int32)
    
# adata <- adata$AnnData(X=X, obs=so@meta.data)
# adata$var_names <- rownames(GetAssayData(so, assay="RNA", slot="counts"))

# adata$raw <- adata
# adata$write_h5ad(h5ad_pp_file)

Session info¶

In [59]:
sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Red Hat Enterprise Linux 8.8 (Ootpa)

Matrix products: default
BLAS/LAPACK: /nobackup/peer/fdeckert/miniconda3/envs/r.4.1.0-FD20200109SPLENO/lib/libopenblasp-r0.3.21.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] viridis_0.6.3         viridisLite_0.4.1     patchwork_1.1.2      
 [4] RColorBrewer_1.1-3    reticulate_1.28       circlize_0.4.15      
 [7] ggplotify_0.1.0       ComplexHeatmap_2.10.0 forcats_1.0.0        
[10] stringr_1.5.0         dplyr_1.1.2           purrr_1.0.1          
[13] readr_2.1.4           tidyr_1.3.0           tibble_3.2.1         
[16] ggplot2_3.4.2         tidyverse_1.3.2       SeuratObject_4.1.3   
[19] Seurat_4.1.1         

loaded via a namespace (and not attached):
  [1] utf8_1.2.3            tidyselect_1.2.0      htmlwidgets_1.6.2    
  [4] Rtsne_0.16            munsell_0.5.0         codetools_0.2-19     
  [7] ica_1.0-3             pbdZMQ_0.3-9          future_1.32.0        
 [10] miniUI_0.1.1.1        withr_2.5.0           spatstat.random_3.2-3
 [13] colorspace_2.1-0      progressr_0.13.0      highr_0.10           
 [16] knitr_1.42            uuid_1.1-0            rstudioapi_0.14      
 [19] stats4_4.1.0          ROCR_1.0-11           tensor_1.5           
 [22] listenv_0.9.0         labeling_0.4.2        repr_1.1.6           
 [25] polyclip_1.10-4       farver_2.1.1          rprojroot_2.0.3      
 [28] parallelly_1.35.0     vctrs_0.6.5           generics_0.1.3       
 [31] xfun_0.39             timechange_0.2.0      R6_2.5.1             
 [34] doParallel_1.0.17     clue_0.3-64           locfit_1.5-9.7       
 [37] spatstat.utils_3.0-4  gridGraphics_0.5-1    promises_1.2.0.1     
 [40] scales_1.2.1          googlesheets4_1.1.0   gtable_0.3.3         
 [43] Cairo_1.6-2           globals_0.16.2        goftest_1.2-3        
 [46] rlang_1.1.1           systemfonts_1.0.5     GlobalOptions_0.1.2  
 [49] splines_4.1.0         lazyeval_0.2.2        gargle_1.4.0         
 [52] spatstat.geom_3.2-9   broom_1.0.4           reshape2_1.4.4       
 [55] abind_1.4-5           modelr_0.1.11         backports_1.4.1      
 [58] httpuv_1.6.9          tools_4.1.0           ellipsis_0.3.2       
 [61] spatstat.core_2.4-4   kableExtra_1.4.0      BiocGenerics_0.40.0  
 [64] ggridges_0.5.4        Rcpp_1.0.10           plyr_1.8.8           
 [67] base64enc_0.1-3       rpart_4.1.23          deldir_1.0-6         
 [70] pbapply_1.7-0         GetoptLong_1.0.5      cowplot_1.1.1        
 [73] S4Vectors_0.32.4      zoo_1.8-12            haven_2.5.2          
 [76] ggrepel_0.9.3         cluster_2.1.4         fs_1.6.2             
 [79] here_1.0.1            magrittr_2.0.3        data.table_1.14.8    
 [82] scattermore_1.0       lmtest_0.9-40         reprex_2.0.2         
 [85] RANN_2.6.1            googledrive_2.1.0     fitdistrplus_1.1-11  
 [88] matrixStats_0.63.0    hms_1.1.3             mime_0.12            
 [91] evaluate_0.20         xtable_1.8-4          readxl_1.4.2         
 [94] IRanges_2.28.0        gridExtra_2.3         shape_1.4.6          
 [97] compiler_4.1.0        KernSmooth_2.23-21    crayon_1.5.2         
[100] htmltools_0.5.5       mgcv_1.8-42           later_1.3.1          
[103] tzdb_0.4.0            lubridate_1.9.2       DBI_1.1.3            
[106] dbplyr_2.3.2          MASS_7.3-58.3         Matrix_1.5-4         
[109] cli_3.6.1             parallel_4.1.0        igraph_1.3.1         
[112] pkgconfig_2.0.3       sp_1.6-0              IRdisplay_1.1        
[115] plotly_4.10.1         spatstat.sparse_3.0-3 xml2_1.3.3           
[118] foreach_1.5.2         svglite_2.1.3         rvest_1.0.3          
[121] yulab.utils_0.0.6     digest_0.6.31         sctransform_0.3.5    
[124] RcppAnnoy_0.0.20      spatstat.data_3.0-4   rmarkdown_2.21       
[127] cellranger_1.1.0      leiden_0.4.3          edgeR_3.36.0         
[130] uwot_0.1.14           shiny_1.7.4           rjson_0.2.21         
[133] lifecycle_1.0.3       nlme_3.1-162          jsonlite_1.8.4       
[136] limma_3.50.3          fansi_1.0.4           pillar_1.9.0         
[139] lattice_0.21-8        fastmap_1.1.1         httr_1.4.5           
[142] survival_3.5-5        glue_1.6.2            png_0.1-8            
[145] iterators_1.0.14      stringi_1.7.6         IRkernel_1.3.2       
[148] irlba_2.3.5.1         future.apply_1.10.0